/**
 * @file Classical_RK.h
 * @brief The class for Classical Runge-Kutta method.
 * @author XDDDD
 * @version 
 * @date 2021-04-28
 */

#ifndef __PAISLEYPARK__CLASSICAL_RK_H__
#define __PAISLEYPARK__CLASSICAL_RK_H__

#include "TimeIntergrator.h"

template <std::size_t DIM>
class Classical_RK : public TimeIntergrator<DIM> {
public:
	Classical_RK(){};
	~Classical_RK(){};
	void one_step_method(std::vector<Point<DIM> >& _Pnts, const double _stepl, VectorFunction<DIM, DIM>& _f) const override;
};

template <std::size_t DIM>
void Classical_RK<DIM>::one_step_method(std::vector<Point<DIM> >& _Pnt, const double _stepl, VectorFunction<DIM, DIM>& _f) const {
	if(_Pnt.size() != 1) {
		std::cerr << "Too many starting points!" << std::endl;
		std::exit(-1);
	}
	
	double time = _Pnt[0].get_time();
	double c[4] = {0, 0.5, 0.5, 1.0};
	double a[4] = {0, 0.5, 0.5, 1.0};
	double b[4] = {1.0/6, 1.0/3, 1.0/3, 1.0/6};
	Point<DIM> y[4];
	for(int i = 0; i != 4; i++) {
		y[i] = _Pnt[0];
		if(i != 0) {
			y[i] = y[i] + a[i]*_stepl*y[i - 1];
			y[i].set_time(y[i - 1].get_time() + c[i]*_stepl);
		}
		y[i] = _f(y[i]);
	}
	for(int i = 0; i != 4; i++) {
		_Pnt[0] = _Pnt[0] + _stepl*b[i]*y[i];
	}
	_Pnt[0].set_time(time + _stepl);
};

#else
//Do nothing.
#endif
